\documentclass[a4paper,10pt]{article}
\usepackage[pdftex]{graphicx}
\usepackage[T1]{fontenc}
\usepackage[slovene]{babel}
\usepackage[utf8x]{inputenc}
\usepackage{hyperref}
\usepackage{amssymb}
\usepackage{amsmath}


\begin{document}
\newcommand{\prm}{r_{M}}
\newcommand{\prp}{r_{P}}
\newcommand{\pqm}{q_{M}}
\newcommand{\pqp}{q_{P}}
\newcommand{\ptt}{\tau}
\newcommand{\pmm}{m}
\newcommand{\pnn}{n}
\newcommand{\pkk}{k}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}

\title{Building a custom genetic oscillator: An alternative way to search the parameter space}
\date{2. 4. 2012}
\author{Martin Stražar}

\maketitle



\section*{Abstract}

\section*{Introduction}

In the following article, we present an alternative method towards building a genetic oscillator. Using a system of two differential equations, we obtain the parameters for a custom oscillator (defined by oscillation frequency and amplitude). That is achieved by the usage of an efficient genetic algorithm and numerical methods.
In \cite{scheper99}, a mathematical model for a single negative feedback loop was proposed. The parameter space was not explored in it's entirety. In this work, we propose an alternative approach to explore the parameter space: given the desired behaviour of the oscillator, the parameters that produce it are calculated.

%TODO: omeni repressilatorje, preveč komplicirano
%TODO: dodaj metode za dosego končnega cilja realizacije
%TODO: poišči primere nadzora velikosti celice, interakcij med proteini in to. Naj se nanaša na konkreten članek z rezultati. 
%TODO: poženi zadevo za kakšen primer v člankih.

%TODO: POTREBEN NADZOR NASLEDNJIH ZADEV:
%0. hitrosti produkcije in degradacije proteinov in mRNA
%1. časovni zamik med transkripcijo in translacijo. (Produciran mRNA in produkcija proteina)
%3. nelinearnost 


%TODO: post-trankripcijska regulacija. SiRNA, nastanejo iz dsRNA z vnosom endonukleaze (Dicer). Potencialni nadzor nelinearnosti v sodelovanju.

%Če je možno dati pod isti promotor:
%P...P...P...(n-krat)-linker-DIC...DIC...DIC...(m-krat)... Nadzor Hillovega koeficienta na večji skali!

%Hillov koeficient: več proteina P -> počasnejša produkcija mRNA.
%Več proteina P -> več endonukleaze -> več siRNA -> manj mRNA. Fuzija proteina P v istem zaporedju (genu) in endonukelaze (ortogonalnost)?

%Nelinearnost vezave: 
%m<0: več mRNA -> več proteina (hitrost narašča s korenom). Logično - v neki točki manj prostih ribosomov?
%m>0: več mRNA -> več proteina (hitrost narašča eksponentno). Popolnoma nelogično? 



%TODO: nadzor hitrosti reakcij v tekočih medijih?

\section*{Methods}\ref{dmdt} and \ref{dpdt}

  \subsection*{Modelling of the oscillator}
    Building the work of \cite{scheper99}, the model of the oscillator is defined as follows: \\
      \beq \frac{dM}{dt} = \frac{r_{M}}{1+(P/k)^n}-q_{M} \cdot M \label{dmdt} \eeq
      \beq \frac{dP}{dt} = r_{P} \cdot M(t-\tau)^m-q_{P} \cdot P \label{dpdt} \eeq


    %TODO: incorporiraj tabelo
\begin{table}[ht!]
    \begin{tabular}{|l|l|l|}
        \hline
        Parameter label & Meaning   & Initial value                                  \\ \hline
        $r_{M}$              & Rate of production (mRNA)                  & 1.0\\ 
        $r_{P}$              & Rate of production (protein)               & 1.0\\
        $q_{M}$              & Rate of degradation (mRNA)                 & 1.0\\ 
        $q_{P}$              & Rate of degradation (protein)              & 1.0\\ 
        $\tau$               & Transcription-translation delay            & 1.0\\ 
        $k$                  & Scaling constant                           & 1.0\\ 
        $m$                  & nonlinearity in protein synthesis cascade  & 3.0\\  
        $n$                  & Hill coeficient                            & 2.0\\
        \hline
    \end{tabular}
\end{table}


    %TODO: Model oscilatorja
    For the purpose of the algorithm, the oscillator is modelled as a static object, defined by an array of parameter values, the running time and the step size. The object model contains and array of points in time and the corresponding values of the observed quantities, calculated using the above differential equations at each step in time.
    Given the above data, it is possible to calculate the amplitude, frequency and wheter the oscillator converges to a steady state or not. These values are used for the purpose of the fitness function, the essential part of the genetic algorithm. 
 

  \subsection*{Initial definitions}  
  %TODO: intitial definitions

  %TODO: Model tarče
  \subsection*{Defining the input}
  The input to the algorithms are the target frequency (in Hz, denoted $\omega$ ) and target amplitude (in units of concentration, denoted $A$). The target function is the represented as a sine wave on the interval of interest:\\
  \beq y = A \cdot sin(\omega x) + \frac{A}{2}. \eeq


  %TODO: opis algoritma
  \subsection*{Finding the parameter values}
  To find the parameter values, that produce a behaviour that best approximate the target sine wave a customied variant of a genetic algorithm is used. A member of the population (an entity) is an oscillator model described above.
  
  \subsubsection*{Evaluating the oscillators}
  %TODO:

  In order to evaluate the behaviuor of the oscillator, the equations \ref{dmdt} and \ref{dpdt} are evaluated using Euler's numerical integration method. During this phase, the sampling interval and the step size should be large- and small enoug respectively to produce accurate enough results. The results are vectors of discrete sample values for concentrations of mRNA and protein in the given sampling interval:
\beq mRNA = [mRNA_{1},...,mRNA_{n}] \label{mrnapoints} \eeq
\beq P = [p_{1},...,p_{n}] \label{ppoints} \eeq

Considering the equations \ref{dmdt} and \ref{dpdt}, a state point of the oscillator in a given time $t$ is then defined as:
\beq s_{t} = [P(t), mRNA(t), mRNA(t- \tau )]\eeq

Along with the computation of \ref{mrnapoints} and \ref{ppoints}, the state points are stored in a vector of states.
\beq S = [s_{i},...,s_{t}] \label{statevector} \eeq

The values of amplitude ($A_{osc}$ - the difference in the peak and low values of \ref{dpdt}) and frequency ($\omega_{osc}$ - the distance from peak to peak) are then derived directly from the values in \ref{dmdt} and \ref{dpdt}. In order to know if a given oscillator is indeed oscillating, the following three condition must be satisfied:
\begin{enumerate}
  \item{$A_{osc} > 0$}
  \item{$\omega_{osc} > 0$}
  \item{There exist two equal state points in \ref{statevector} and the same state point value is not repeated for a longer interval than $\tau$}
\end{enumerate}

  \subsubsection*{The fitness evaluation function}
  The quality of an entity is measured using the following metrics:
%  \begin{itemize}
%    \item{The amplitude of the oscillations,}
%    \item{The frequency of the oscillations,}
%    \item{The symmetry in rise and fall times.}
%  \end{itemize}

\begin{table}[ht!]
    \begin{tabular}{|l|l|l|}
        \hline
        Shorthand & Metric & Evaluation \\ \hline
        $M_{A}$ & $Amplitude$                   & $|A_{sin} - A_{osc}|$                               \\ 
        $M_{\omega}$ & $Frequency$              & $|\omega_{sin} - \omega_{osc}|$                               \\ 
        $M_{sym}$ & $Symmetry$                  & $|t_{rise} - t_{fall}|$          \\ 
        \hline
    \end{tabular}
\end{table}

  For all the mentioned metrics, the smallest error in comparison with the target sine wave is considered best. Summing up the achieved values, the target fitness function is defined as the mean of the subgroup of the above metrics:\\
\beq F_{osc} = \frac{M_{A} + M_{\omega} + M_{sym}}{3}\eeq
Hence, the lower the value of the fitness function, the better the approximation to the target sine wave.

\subsubsection*{Population dynamics}
%TODO: opiši dinamiko populacije
The algorithm is run with a population of oscillators with a starting size. In each generation, the parameters are modified and the best individuals are selected, depending on their $F_{osc}$ values. According to the mean fitness value in the population, the oscillating individuals with the value lower than the mean progress in to the next generation and produce offspring. The resulting size of the remaining population is then at most half the size of the original population:
\beq F_{mean} = E(\sum_{osc} F_{osc}) \eeq
\beq Pop_{rem} = |\{osc;  F_{osc} <= F_{mean}\}| \eeq

\subsubsection*{Offspring production}
A child oscillator $Osc_{c}$ is an oscillator, derived from the original oscillator $Osc_{p}$. During the derivation process, the predetermined value for a mutation %M% is used.
\beq mut = P_{mut} \cdot rnd_{1} \cdot rnd_{2}, \eeq
where $rnd_{x}$ denotes a random float value in range $0$ to $1$ and $P_{mut}$ the predetermined probability of an occuring mutation. A child oscillator is then defined as:\\
\beq Osc_{p}:[\ptt,\prm,\prp,\pqm,\pqp,\pmm,\pnn,\pkk] \rightarrow \eeq
\beq Osc_{c}:[\ptt+mut,\prm+mut,\prp+mut,\pqm+mut,\pqp+mut,\pmm+mut,\pnn+mut,\pkk+mut] \eeq \\\\

The offspring is defined as a number of child oscillators, derived from an original parent oscillator in each generation. Parent oscillators with a lower fitness value get to have a larger part of the total offspring:
\beq Off_{total} = 1 + Pop_{rem} \eeq
\beq Off_{osc} = \frac{1}{F_{osc}} \cdot \frac{1}{\sum_{i=0}^Pop_{rem} F_{i}} \cdot Off_{total}. \eeq




  \subsubsection*{The generation loop}
    The starting point of the algorithm is an array of parameters (given in table 1), which determines the initial population. By adjusting the values of parameter array, new members of population are born. A fitness function is used in each generation iteration to determine the survival of the fittest, by sorting the population from best to worst match. Depending on each individuals' fitness value, the number of descendants is calculated and the offspring is produced. Consequently, the best match for the sine wave if found.
    The generation loop is repetead for a predetermined (experimentally derived) number of iterations. 

  The generation loop is repeated in three steps, where each of the steps uses a different subgroup of the metrics to achieve its subgoal. The output (achieved parameters) of each step is the input of the step following it. Table 2 describes the metrics that are evaluated in each of the three steps.

\begin{table}[ht!]
    \begin{tabular}{|l|l|l|}
        \hline
        Step & Evaluated metrics & $P_{mut}$ \\ \hline
        $0$              & Amplitude & $\frac{1}{8}$                              \\ 
        $1$              & Frequency & $\frac{1}{8}$                              \\ 
        $2$              & Amplitude, Frequency, Symmetry & $\frac{1}{8}$          \\ 
        \hline
    \end{tabular}
\end{table}

 
S




\section*{Test results}


\section*{Discussion}
%TODO: subject to finding the proper parameter values
%TODO: by finding the relevant sources,  the correct paramter values could be adjusted

\section*{References}
\begin{thebibliography}{9}
\bibitem{scheper99}
   Tjeerd olde Scheper, Don Klinkenberg, Cyriel Pennartz and Jaap van Pelt
   \emph{A Mathematical Model for the Intracellular Circadian Rhythm Generator} 
   The Journal of Neuroscience,
   1 January 1999, 
   19(1):40-47

\bibitem{man99}
  K. F. Man, K. S. Tang and S. Kwong,
  \emph{Genetic algorothms: concepts and design},
  Springer, 
  1999, 
	1-85233-072-4

 
\end{thebibliography}

\end{document}

